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SUMMARY 


This  report  presents  an  iterative  method  for  inverting  VLF/LF 
ionosound  data  to  obtain  ionospheric  conductivity  profiles.  The  re¬ 
lationship  between  the  detail  that  can  be  realized  in  the  calculated 
profile  and  the  quality  and  quantity  of  the  reflection  data  is  given 
explicitly  by  our  approach.  We  consider  only  altitudes  below  about 
70  km  where  the  propagation  can  be  assumed  isotropic,  although  the 
method  can  be  extended  to  anisotropic  propagation. 

We  demonstrate  the  method  by  applying  it  to  two  qualitatively 
different  ionospheric  models — a  strong  SPE  disturbance  and  a  C-layer 
exhibiting  a  well-defined  conductivity  peak.  Ground-level  reflection 
data  contain  information  about  only  those  heights  from  which  signifi¬ 
cant  energy  is  returned.  Our  calculated  profiles  agree  well  with 
the  true  profiles  over  the  height  range  where  the  strongest  reflec¬ 
tions  occur,  and  break  down  outside  of  this  range.  The  solutions 
converge  toward  the  true  profiles  despite  our  use  of  initial  esti¬ 
mates  that  contain  no  information  about  the  true  ionospheric  struc¬ 
ture.  The  height  resolution  actually  achieved  for  the  two  examples 
is  much  better  than  the  theoretically  predicted  limit. 
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PREFACE 


This  report  describes  a  method  for  inverting  VLF/LF  reflection 
data  to  obtain  ionospheric  conductivity  profiles  under  conditions 
where  the  propagation  can  be  assumed  isotropic.  Me  demonstrate  the 
method  by  inverting  computer-generated  reflection  coefficients  for 
two  generic  model  ionospheres.  In  addition,  we  present  a  means  of 
choosing  the  best  compromise  between  the  deleterious  effects  of  (1) 
noisy  data  and  (2)  loss  of  uniqueness  due  to  incomplete  data. 
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I.  INTRODUCTION 


Radio-sounding  at  VLF  (3  to  30  kHz)  and  LF  (30  to  300  kHz)  is 
often  used  to  sample  the  ionosphere  at  heights  below  about  100  km 
where  high-frequency  (HF)  ionosounds  provide  little  information. 

Further,  at  altitudes  below  about  70  km  the  electron  density  is  often 
too  small  to  be  sensed  by  rocket-borne  probes,  and  VLF/LF  sounding 
is  the  only  sampling  means  available.  Once  ground-level  VLF/LF  re¬ 
flection  coefficients  have  been  measured,  they  must  be  inverted  to 
infer  ionospheric  structure.  This  "inverse"  problem  is  much  more 
difficult  than  the  "forward"  problem  of  calculating  reflection  coef¬ 
ficients  from  a  specified  model  ionosphere. 

Inversion  of  ionosound  data  is  often  performed  by  trial  and 
error,  with  reflection  coefficients  being  calculated  from  a  number  of 
model  ionospheres,  and  the  model  giving  the  coefficients  that — in  the 
judgment  of  the  analyst — most  closely  agree  with  experiment  being  deemed 
correct.  This  report  describes  a  method  tor  inverting  VLF/LF  reflec¬ 
tion  data  to  obtain  ionospheric  conductivity  profiles  without  using 
trial  and  error.  We  consider  only  altitudes  below  about  70  km  where 
the  propagation  can  be  assumed  isotropic,  although  the  method  can  be 
extended  to  anisotropic  propagation.  This  restriction  on  altitude 
limits  the  applicability  to  disturbed  conditions — such  as  solar- 
proton  events  (SPE) — and  some  ambient  daytime  situations. 

Our  procedure  extends  the  spectral  method  described  by  Parker 
[1977]  for  implementing  the  Backus  and  Gilbert  theory  [1967,  1968, 

1970,  1971],  which  was  developed  for  inverting  seismic  or  electro¬ 
magnetic  soundings  of  the  earth's  interior.  We  find  that  method  to 
be  well  suited  to  the  ionosphere  as  well.  Chuang  and  Yeh  [1977] 
l  applied  the  Backus  and  Gilbert  theory  to  HF  ionosound  data.  However, 

they  used  geometric  optics,  whereas  VLF/LF  propagation  must  usually 
be  described  by  full-wave  theory.  Shellman  [1973]  and  Morfitt  and 
Shellman  [1977]  inverted  VLF/LF  ionospheric  reflection  data,  including 
geomagnetic  anisotropy,  by  performing  a  least-squares  optimization  of 
the  electron  density  and  collision  frequency.  That  approach,  although 
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direct,  requires  that  the  number  and  location  of  the  sampling  heights 
be  selected  beforehand.  The  Backus-Cilbert  method  requires  no  such 
preselection  and  imposes  fewer  constraints  on  the  functional  form  of 
the  inverse  solution.  Moreover,  it  gives  the  relationship  between  the 
detail  that  can  be  realized  in  the  calculated  conductivity  height 
profile  and  the  quality  and  quantity  of  the  reflection  data. 

Section  II  describes  the  mathematical  basis  of  our  inversion 
method.  Section  III  applies  the  method  to  actual  and  computer- 
synthesized  reflection-coefficient  data  and  evaluates  uncertainties 
in  the  inferred  profiles  caused  by  noise  and  incomplete  data.  Sec¬ 
tion  III  also  presents  a  means  of  estimating  the  height  range  over 
which  waves  at  the  sounding  frequencies  interact  strongly  with  the 
ionosphere;  and,  therefore,  the  height  range  over  which  the  sounding 
data  can  provide  information  on  ionospheric  structure.  Section  IV 
gives  the  conclusions,  and  the  Appendix  presents  certain  mathemati¬ 
cal  details. 
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II.  INVERSE  THEORY  APPLIED  TO  IONOSOUNDING  DATA 

Inferring  ionospheric  height  profiles  from  ground-based  reflec¬ 
tion  data  is  one  of  a  class  of  problems  treated  by  inverse  theory,  the 
object  of  which  is  to  use  limited  and  usually  noisy  data  to  determine 
the  maximum  information  about  geophysical  structure.  Inverse  problems 
are  typically  much  more  difficult  than  forward  problems,  which  require 
calculation  of  the  signal  emerging  from  specified  propagation  media. 

This  difficulty  is  due  primarily  to  the  nonuniqueness  of  the  inferred 
profile.  Backus  and  Gilbert  [19671  showed  that — under  fairly  general 
conditions — a  finite  data  set  is  consistent  with  an  infinite  number 
of  possible  profiles,  if  it  is  consistent  with  any  one.  In  addition, 
in  many  cases  (including  the  one  addressed  here),  the  measured  signal 
has  a  nonlinear  dependence  on  the  propagation  medium.  This  nonlinearity 
could  introduce  nonuniqueness,  because  the  unknown  profile  must  be 
constructed  iteratively  from  a  starting  estimate  and  different  start¬ 
ing  estimates  can,  in  principle,  produce  different  profiles. 

Despite  these  difficulties,  the  importance  of  obtaining  informa¬ 
tion  about  the  earth  or  its  environment  from  remote  probes  has  led  to 
at  least  a  partial  resolution  of  questions  pertaining  to  uniqueness 
[Backus  and  Gilbert,  1967,  1968,  1970;  Oldenburg,  1978).  The  non¬ 
uniqueness  due  to  incomplete  data  is  fairly  well  understood  and  char¬ 
acterized.  That  due  to  nonlinear  dependence  on  the  medium  is  still 
unresolved,  but  seems  to  cause  no  difficulty  for  the  numerical  examples 
that  we  give  in  Sec.  III. 

FRECHET  KERNELS  FOR  ISOTROPIC  IONOSPHERIC  CONDUCTIVITY 

The  first  step  in  the  inverse  solution  is  to  characterize  the 
dependence  of  the  reflection  coefficients  on  the  structure  of  the 
ionosphere.  For  TM  (vertically  polarized)  waves  the  reflection  coef¬ 
ficient  R  below  the  ionosphere  can  be  defined  as 


R  =  H(d)/H(u) 

y  /  y 


(i) 
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where  and  are  incident  and  reflected  magnetic  intensities 
measured  at  the  same  point,  and  y  is  the  cartesian  coordinate  normal 
to  the  plane  of  incidence.  The  equation  governing  R  for  plane  waves 
reflected  from  an  isotropic  ionosphere  is  [Budden,  1961] 


2i  dR  .  2..  ,2 

IdI=Cn(1-R) 


-V  (1  +  R)2  , 

n  C 


(2) 


2  2  2 

where  n(z)  is  the  complex  refractive  index,  q  =  n  -  S  ,  S  and  C 
are  the  sine  and  cosine  of  the  incident  wave  in  free  space,  and  k  is 
the  free-space  wave  number.  An  analogous  equation  exists  for  the  TE 
(horizontally  polarized)  reflection  coefficient  [Budden,  1961],  but 
we  will  consider  only  vertically  polarized  waves. 

The  general  refractive  index  expression,  including  geomagnetic 
effects,  is 


2  _  .  e 
n  =1 - 

uie. 


N  (z) 

JL 


in^zMw  -  iV^(z)  +  o)c^(z)] 


(3) 


where  the  sum  is  over  all  charged  species  (electrons,  positive  and 
negative  ions),  a)  is  the  angular  frequency,  e  is  the  electronic  charge, 
Eg  is  the  permittivity  of  free  space,  and  N^,,  m^,  v^,  and  are  the 
height-dependent  number  density,  mass,  collision  frequency,  and  gyro- 
frequency  of  the  ylh  species.  The  isotropic  equation  (2)  is  valid 
provided  all  important  reflections  occur  below  altitudes  of  70  to 
75  km,  where  V  >  u>  .  Moreover,  since  u>  <<  w  <  V  at  VLF/LF  in  this 

Vy  2  w 

low-altitude  regime,  n  assumes  the  simplified  form 


n2 (z)  =-  1  -  i  , 

WE0 


(4) 


where  the  total  conductivity  is  given  by 


0(z)  =-  e 


m  V 
Y  Y 


(5) 
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Heavy  ions  make  important  contributions  to  the  conductivity  only  at 
altitudes  below  50  to  55  km  above;  which  the  conductivity  may  be  con¬ 
sidered  entirely  due  to  electrons.  At  altitudes  where  Eqs.  (2)  and 
(4)  are  valid,  reflection  data  can  provide  information  only  about  the 
total  conduct ivity .  Auxiliary  information  is  needed  to  separate  the 

components  appearing  in  Eq.  (5). 

2 

Equation  (2) — with  n  given  by  Kq .  (4) — can  be  integrated  numer¬ 
ically  for  a  given  u)  to  obtain  R(0)  at  the  ground.  We  start  the  inte¬ 
gration  from  the  initial  estimate 


R(z  > 

m 


n2C  -  q 
n20  +  q 


(t>) 


where  z  is  chosen  large  enough  for  R  to  be  slowlv  varying,  and  Eq . 
m 

(6)  results  from  setting  dR/dz  =  0  in  Eq.  (2).  This  solution  for  the 
reflection  coefficient  R(0) ,  corresponding  to  a  given  conductivity 
profiLe  o( z),  comprises  the  forward  problem. 

The  inverse  problem  is  to  obtain  the  best  possible  estimate  of 
a(z)  from  R(0) ,  specified  at  a  finite  set  of  frequencies.  We  con¬ 
struct  O  by  computing  successive  corrections  to  an  initial  estimate 
0^(z),  each  correction  being  determined  by  relating  a  small  change 
So(z)  in  conductivity  at  height  z  to  a  change  SR^fO)  in  the  ground- 
level  reflection  coefficient  at  frequency  f.  =  <jj./27T.  We  find  the 

l  i 

necessary  relationship  by  expanding  an  equation  analogous  to  Eq .  (2) 
for  R'(z)  =  R(z)  +  <5r(z)  in  terms  of  a’( z)  =  0(z)  +  60(z),  retaining 
only  terms  linear  in  5R  and  5o,  and  using  Eq.  (2).  The  resulting 
linear  first-order  differential  equation  for  6R(z)  is 


_d_ 

dz 


6R  +  P(z)6R(z) 


Q(z)  , 


(7) 


where 


k  r 

i 


Cn  (1  -  R)  + 


A 

n2C 


(1  +  R) 


P(z) 


(8) 
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1  fvo  .  r)2  _  jL  (1  +  M2  to(z)  • 

«*)  •  -  2  ^  [C(1  '  A  J 


Equation 


(7)  can  be  integrated  to  give 


6R(0)  = 


z 

I 

Z 

1  dz  Q(z)  exp 

if  dz' 

P(z’) 

J 

0 

[/o 

“ 

»  2 

1 

5R(z  )  exp 

fll 

[1 

dz  p(z)  ] 

real  component; 

so  the  second  t« 

(9) 


(10) 


it  the  ground,  and 


6R(U) 


dz  G(z)6a(z)  , 


(ID 


where  C  is  the  i-rechet 
ductivity  at  altitude 
ground : 


kernel,  which  relates  a  small  change  in  con 
z  to  a  change  in  reflection  coefficient  at  the 


(z) 


(1  +  R) 


■1 


x  exp 


(1  +  R) 


(12) 
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CONSTRUCTION  OF  THE  INVERSE  SOLUTION 

As  emphasized  by  Backus  and  Gilbert  [1967,  1968,  1970],  the  esti¬ 
mation  of  an  unknown  propagation  medium  from  a  finite  set  of  field 
measurements  requires  two  steps.  First,  we  must  construct  a  conduc¬ 
tivity  profile  that  reproduces  the  data.  Second,  and  tqually  important, 
we  must  characterize  the  class  of  solutions  to  which  our  profile  and — 
presumably — the  true  profile  belong.  Here  we  address  the  first  step. 

The  uniqueness  of  our  inverse  solution  is  considered  in  the  following 
subsect  ion. 

To  construct  the  inverse  solution  we  must  find  a  conductivity 
profile  that  matches  a  set  of  reflection  coefficient  data  at  a  number 
(n)  of  frequencies.  Following  Oldenburg  [1978],  we  pick  an  initial 
estimate  compute  R^*  from  Eq .  (2)  at  each  frequency  f^,  and  set 

SR*  =  R^ata  -  R*?(0)  ,  i  =  1,  ....  n  ,  (13) 

where  superscripts  denote  the  order  of  the  correction  (i.e.,  the 
number  of  iterations),  and  subscripts  index  the  frequencies.  We  then 
express  the  first-order  correction  to  O  ,  6a  ,  as 


6a 


*(z)  =  ^  ^  a.G*(z)  , 

j=l 


(14) 


where  the  ou  are  found  by  solving  the  n  *  n  set  of  equations 


r.,o.  =  6r: 

ij  J  1 


(15) 


i  =  l 


and 


/m 

dz  G1(z)G*(z) 


(16) 
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Next,  K.  is  computed  from  Eq.  (2)  using  the  new  conductivity  estimate 
1  1  0  1 

O  (z)  =  0  (z)  +  60  (z),  and  the  process  is  successively  repeated.  At 
each  iteration  N  we  use  the  rms  relative  error 

X  i=l 


as  a  measure  of  convergence. 

In  practice,  we  use  tne  following  modifications  to  the  above 
procedure : 

1.  Because  the  data  and  Frechet  kernels  are  expressed  in  com¬ 
plex  form,  the  differential  conductivities  So1  produced  by 
Eq.  (14)  can  have  imaginary  parts.  We  have  found  it  suffi¬ 
cient  to  use  the  fact  that  the  conductivity  in  Eq.  (14)  must 
be  real  and  simply  suppress  the  imaginary  part  of  the  cor¬ 
rection  at  each  iteration.  A  more  rigorous — and  complicated — 
approach  would  be  to  express  all  quantities  in  real  form  at 
the  outset. 

2.  As  discussed  by  Oldenburg,  direct  use  of  the  output  correc¬ 
tions  to  the  conductivity  at  each  iteration  tends  to  pro¬ 
duce  oscillatory  profiles,  and  a  nonphysical  build-up  of 
conductivity  at  low  heights.  These  pathologies  require  that 
each  correction  be  substantially  smoothed  and  filtered  before 
being  used  in  the  next  iteration.  In  addition,  we  do  not 
allow  the  total  conductivity  at  any  iteration  to  become 
negative  but  redefine  it  as  zero,  if  required,  before  pro¬ 
ceeding.  These  steps  slow  the  convergence,  but  appear 
necessary  to  produce  realistic  profiles. 

3.  Following  Parker  [  1 9 7 7  |  we  diagonalize  the  matrix  F  before 
inversion,  so  the  expansion  | Eq .  (14)]  is  actually  performed 
in  the  eigenfunctions  of  the  T  matrix.  This  step  allows  each 
term  in  the  diagonalized  basis  to  be  associated  with  a 
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relative  uncertainty  in  the  data,  and  components  that  intro¬ 
duce  large  statistical  uncertainties  can  he  omitted.  The 
Appendix  gives  the  details  of  this  spectral  expansion  method. 

UNIQUENESS  OF  THE  INVERSE  SOLUTION 

The  above  algorithm  is  iterated  until  the  relaf;ve  rms  error 
[Eq.  (17)]  asymptotically  approaches  a  minimum.  We  have  then  done  as 
well  as  possible  in  constructing  a  conductivity  profile  with  the  data 
and  starting  estimate  used.  The  remaining  task  is  to  assess  the 
validity  of  the  profile  at  various  heights.  As  discussed  above  there 
are  three  sources  of  nonuniqueness: 

1.  Noise  in  the  reflection  coefficient  data  causes  uncertainties. 

2.  The  data  set  is  finite,  and  therefore  incomplete. 

3.  The  conductivity  profile  and  reflection  coefficient  data  are 
nonlinearly  related. 

The  first  error  source  is  easily  assessed  with  multivariate 
error  analysis.  The  spectral  expansion  method  (see  Appendix)  provides 
a  simple  means  of  estimating  the  effect  of  measurement  uncertainties 
on  the  accuracy  of  the  inferred  conductivity  profile.  More  important, 
variances  in  the  reflection  data  are  included  in  the  solution  in  a 
manner  that  deemphasizes  frequencies  at  which  measurements  are  ex¬ 
ceptionally  noisy. 

The  second  type  of  nonuniqueness  has  received  the  most  attention 
in  the  literature,  since  it  is  present  regardless  of  the  care  taken 
to  reduce  experimental  uncertainty.  It  is  therefore  the  ultimate 
limitation  on  using  a  finite  data  set  to  characterize  an  unknown  func¬ 
tion.  Before  specifically  addressing  this  type  of  nonuniqueness,  we 
consider  a  somewhat  similar  situation,  which  is  helpful  in  understand¬ 
ing  inaccuracies  caused  by  incomplete  data. 

The  Fourier  series  representation  of  an  arbitrary  periodic  func¬ 
tion  <f>(x) ,  0  <  x  <  2v ,  is  given  by 
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<1>(x) 


1 

/2-n 


<JO 


c 

n 


inx 

e 


(18a) 


where 


C 

n 


dx  <{>(x) 


-inx 

e 


(18b) 


It  is  well  known  that — subject  to  some  weak  requirements — the  infinite 
set  of  coefficients  {C  }  completely  defines  <p.  If,  however,  only  a 
finite  set  -N  £  n  £  N  of  Fourier  coefficients  is  given,  the  high- 
frequency  content  of  <p  will  be  lost;  and  any  <£'  that  differs  from  (f 
only  in  these  neglected  components  will  have  the  same  coefficients 
C  ,  for  ] n I  £  N. 

To  quantify  the  uncertainty  caused  by  omitting  the  coefficients 
of  index  |n|  >  N,  we  let  (<t>)  denote  the  estimate  for  <J>  based  on  the 
incomplete  set  of  C  ,  and  write 

<+(*»  -  -f  V  Cn  e‘“  .  (19) 

SR 


By  using  Eq.  (18b),  we  find 


<<J>(x)) 


$(x')An(x,  x')  , 


(20a) 


where 


yx,  x')  =  ±  ^ 


in(x-x' ) 


n=*-N 


(20b) 
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Ajj  is  the  Dirichlet  kernel,  and  is  a  measure  of  how  closely  <p  can  be 
resolved  by  the  finite  harmonic  set.  The  width  of  A^  is  proportional 
to  1/N,  so  the  resolution  improves  as  the  number  of  "data"  points 
becomes  larger.  In  the  limit  N  -*•  00 ,  (<j>)  -*■  <p,  and  the  resolution  be¬ 
comes  perfect.  This  behavior  occurs  because  A^  approaches  a  Dirac 
delta  function  for  large  N. 

The  above  approach  can  be  used  in  the  ionospheric  inverse  problem 
to  quantify  the  loss  of  resolution  due  to  having  incomplete  reflection 
data.  For  a  given  iteration,  Eq.  (14)  can  be  written 


(6a(z)>  =  ^  a^G* (2)  , 

j=l 


(21) 


where  we  place  (  )  around  60  to  emphasize  that  only  a  finite  set  of 
frequencies  is  available,  and  our  estimate  can  therefore  differ  from 
the  "true"  correction  6a.  By  using  Eqs.  (15)  and  (11),  we  find 


a 


J 


n 


k*l 


«Rk(0) 


dz'  Gk(zT)6o(z’) 


(22) 


I 


and  substitution  into  Eq.  (21)  gives 


<«o(z)) 


dz'  6o(z')A(z,  z')  , 


(23) 


with 
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n 

A(Z,  z’)  =  rjk  Cj(z)Gk(z,)  •  (24> 

j,k=l 

The  averaging  function  [Eq.  (24)]  is  the  analog  of  the  Dirichlet  kernel 
for  Fourier  series.  Variations  of  the  conductivity  on  a  scale  shorter 
than  the  width  of  A(z,  z')  about  a  given  height  will  not  be  resolvable 
by  the  data  set  used  to  generate  the  profile.  Even  though  we  con¬ 
structed  the  averaging  function  for  a  small  variation  in  conductivity, 

2 

Backus  and  Gilbert  [1968]  prove  that — to  order  | 6o  (  z )  |  — Eq .  (24)  pro¬ 
vide-,*  i  measure  of  the  height  resolution  of  the  total  conductivity. 

show  in  the  Appendix  that  Eq .  (24)  can  be  rewritten  in  essen¬ 
tially  ..'ic  form 


k«sn 

A(z ,  z’)  =  H*(z)H.(z')  ,  (25) 

i=l 


where  the  functions  H ^ (z )  are  produced  by  orthonormalizing  the  set  of 
Frechet  kernels  G^(z).  As  for  the  orthonormal  Fourier  series  basis 
function  discussed  above,  the  H.  functions  can  be  ranked  in  order  of 

l 

increasingly  oscillatory  behavior.  The  diagonal  form  of  the  averag¬ 
ing  function  given  by  F.q .  (25)  is  convenient  for  evaluating  the  im¬ 
portance  of  retaining  the  more  oscillatory  components.  In  general, 
a  better  fit  to  the  data  is  achieved  by  increasing  the  number  of  com¬ 
ponents  used  in  the  construction  of  <Sa(z).  However,  the  resulting 
conductivity  profile  can  exhibit  nonphysical  oscillations  if  too  many 

-i. 

terms  are  retained.  The  optimum  number  of  components  depends 
partly  on  the  estimated  noise  in  the  data,  and  partly  on  the  analyst's 
intuition  regarding  profile  smoothness.  The  relation  of  resolution 

This  situation  is  analogous  to  that  which  can  be  encountered  in 
fitting  a  high-order  polynomial  to  a  large  number  of  data  points.  The 
polynomial  can  match  the  data  exactly  at  each  point,  but  exhibit 
spurious  oscillations  between  the  points. 
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to  noise  is  treated  in  the  next  section  and  the  Appendix.  Its  rela¬ 
tion  to  the  analyst's  intuition  is  difficult  to  quantify,  however. 

The  remaining  source  of  nonuniqueness  is  the  nonlinear  relation 
between  conductivity  and  the  ground-level  reflection  coefficient. 

This  source  is  by  far  the  most  difficult  to  characterize  in  general, 
because  there  is  no  guarantee  that  other  "globally"  different  conduc¬ 
tivity  profiles  would  not  be  obtained  with  different  initial  esti¬ 
mates  cr°.  The  most  practical  means  of  exploting  this  possibility  is 
to  choose  widely  different  starting  estimates  and  see  if  the  calcu¬ 
lated  profiles  agree,  aside  from  minor  variations.  We  have  found  that 
our  algorithm  produces  iterative  solutions  even  with  a  null  initial 
estimate  (0°  =  0).  Moreover,  calculations  made  from  the  same  data 
sets — but  with  various  nonzero  estimates — gave  virtually  identical 
profiles.  These  results  strongly  suggest  (but  do  not  prove)  that — 
in  practice — the  nonlinear  dependence  of  the  reflection  coefficients 
on  ionospheric  structure  causes  no  loss  of  uniqueness. 
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III.  NUMERICAL  EXAMPLES 


This  section  applies  the  inversion  method  developed  in  Sec.  II 
and  the  Appendix.  Recall  that  our  assumption  of  isotropic  propagation 
restricts  this  application  to  conditions  where  all  important  reflec¬ 
tions  occur  below  70  to  75  km.  We  use  the  examples  given  below  to 
test  the  ability  of  our  inversion  method  to  obtain  correct  conduc¬ 
tivity  height  profiles  from  a  set  of  ground-level  reflection  coef¬ 
ficients;  and  to  assess  the  effects  of  measurement  noise  and  incom¬ 
plete  data  on  the  accuracy  and  height  resolution  of  our  results.  The 
examples  were  chosen  to  illustrate  the  salient  aspects  of  our  approach 
and  are  intentionally  simple.  Our  method  should  work  equally  well  for 
more  complicated  situations. 

STRONG  SPE 

Our  first  example  is  based  on  high-latitude  satellite  measure¬ 
ments  of  incident  particle  fluxes  at  1508  UT  time  on  4  August  1972. 
These  measurements  were  made  near  the  peak  of  one  of  the  strongest 
SPF.s  on  record.  Height  profiles  of  electron  and  ion  densities  were 
calculated  from  the  measured  fluxes  with  Lockheed's  air-chemistry  code, 
and  provided  to  the  authors  by  Reagan  [1978].  These  densities,  along 
with  nominal  height  profiles  of  electron  and  ion  collision  frequencies 
and  masses  were  used  with  Eq.  (5)  to  obtain  the  conductivity  height 
profile  shown  (solid  line)  in  Fig.  1.  This  profile  is  an  example  of 
an  extremely  severe  natural  disturbance. 

Here  we  illustrate  application  of  our  method  rather  than  perform 
detailed  calculations  for  specific  situations.  Therefore,  for  conve¬ 
nience,  we  approximate  the  strong  SPE  profile  with  the  often-used  form 


where  *  3.94  >  10  J ^  mhos/m  and  0  =  0.39  km  Figure  1  shows  the 
agreement  between  the  approximation  [Eq.  (26)]  and  the  actual  conduc¬ 
tivity  profile. 


Inverse  Solution 


To  provide  a  starting  point  for  our  inverse  solution,  we  have 
calculated  TM  reflection  coefficients  for  the  profile  given  by  Eq. 

(26)  by  integrating  Eq.  (2)  for  an  incidence  angle  of  65  deg  at  fre¬ 
quencies  of  4,  6,  9,  14,  21,  30,  and  48  kHz.  Our  reason  for  choosing 
this  set  of  frequencies  is  discussed  below.  These  coefficients  serve 
as  artificial  "data,*’  which  are  analogous  to  measured  data  that  would 
be  available  from  ionosounders .  To  suppress  the  unphysical  build-up 
of  conductivity  below  35  km,  we  use  the  weighting  technique  of  Olden¬ 
burg  [1978]  (see  Appendix),  and  choose 

(  0  ,  z  £  b 

v(z)  =  <  (27) 

(  1  -  exp  |-Y(z  -  b) J  ,  z  ^  b 

for  the  weighting  function  w(z)  with  b  =  35  km,  y  =  0.5  km  1 .  We 
compute  successive  corrections  to  the  conductivity  using  the  spectral 
expansion  method  and  a  null  (a®  =  0)  starting  estimate.  Therefore, 
our  initial  profile  contains  no  information  about  the  true  profile. 
Oscillations  in  the  corrections  6o  for  each  iteration  are  controlled 
by  our  use  of  a  weighted  average  of  the  correction  at  height  z  (weight 
0.5)  and  at  0.7  km  above  and  below  z  (weight  0.25). 

Figure  2  shows  the  computed  conductivity  profile  following  the 
1,;,'  ,  5th,  10!!.,  and  1 5th  iterations,  the  true  profile  (dashed  lines), 
and  the  relative  rms  error  e.  The  computed  profile  was  still  slowly 
increasing  near  50  km  after  the  15 th  iteration,  and  more  iterations 
presumably  would  have  given  closer  agreement  with  the  true  profile 
between  49  and  50  km.  Figure  3  gives  an  expanded  plot  of  the  calcu¬ 
lated  and  true  profiles  after  the  15 th  iteration.  The  agreement  is 
good  at  altitudes  between  about  40  and  49  km.  Outside  of  this  alti¬ 
tude  range  the  calculated  conductivity  is  much  smaller  than  the  true 
conductivity.  We  give  the  physical  reason  for  this  disagreement  below. 

Sen sitivity  to  Noise  and  Limited  Data 

To  evaluate  the  effect  of  measurement  noise  on  the  inferred  con¬ 
ductivity,  we  include  hypothetical  uncertainties  Sj  in  the  reflection 


Conductivity  (mhos/m) 
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Height  (km) 

Fig.  3--Calculated  and  true  profiles  after  15  iterations  for  the  strong  SPE 
example.  Vertical  bars  show  uncertainty  caused  by  hypothetical  noise 
in  the  reflection  coefficients.  Horizontal  bars  show  uncertainty 
caused  by  incompleteness  of  the  reflection  data. 
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coefficients  used  in  our  calculations.  The  values  for  are  derived 
from  typical  measurement  uncertainties  (Rasmussen  [1981])  and  are 
listed  in  Table  1. 


Table  1 

ASSUMED  UNCERTAINTIES  IN  REFLECTION  COEFFICIENTS 
FOR  THE  STRONG  SPE  EXAMPLE 


Frequency 

(kHz) 

!s./r.(0) | 

4 

0.3 

b 

0.3 

9 

0.3 

14 

0.25 

21 

0.1 

30 

0.1 

48 

0.15 

By  including  S^,  we  can  examine  sensitivity  of  the  profile  at  various 
heights  to  me;  urement  noise,  even  though  the  computer-generated  "data" 
are  noiseless. 

The  variance  in  conductivity  is  defined  as 

Var  (a(z) ]  =  <a2(z)>  -  <a(z))2  ,  (28) 

where  (  )  denotes  averages.  The  Appendix  shows  that 

kgn 

Var  [o(z)]  =  fw(z)]2y^X~  |H.(z)|2  ,  (29) 

i=l  i 

where  the  sum  is  over  the  positive-definite  contributions  from  each 
diagonalized  component  H^(z).  The  eigenvalues  X^  for  the  15 th  iter¬ 
ation  are  listed  in  Table  2,  as  is  the  estimated  variance  at  height 
44  km  formed  by  truncating  F.q.  (29)  at  different  values  of  k.  Table 
2  shows  that  retention  of  fewer  terms  in  Eq.  (29)  causes  a  decreased 


T  - 
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Table  2 

EFFECTS  OF  NOISE  FOR  STRONG  SPE  EXAMPLE 


Cumulative 


[Var  (a)]172 

Eigenvalue  (mhos/m) 


i 

— 

1 

1.81 

X 

10~6 

k 

3* 

1 

1.70 

X 

10~9 

i 

= 

2 

4.32 

X 

Iff  7 

k 

= 

2 

3.89 

X 

io'9 

i 

= 

3 

9.46 

X 

10' 9 

k 

= 

3 

1.69 

X 

00 

1 

o 

i 

= 

4 

1.72 

X 

10~9 

k 

= 

4 

2.38 

X 

10  8 

i 

= 

5 

4.51 

X 

o 

rH 

1 

o 

k 

ss 

5 

9.19 

X 

00 

1 

C 

•— i 

i 

= 

6 

1.95 

X 

10'11 

k 

= 

6 

2.35 

X 

10  7 

i 

= 

7 

2.65 

X 

io~13 

k 

= 

7 

5.40 

X 

io'6 

variance.  If  variance  were  the  only  consideration,  retention  of  only 
the  first  spectral  component  would  provide  the  "best"  solution.  How¬ 
ever,  as  discussed  below,  there  is  a  tradeoff  between  noise  variance 
and  height  resolution,  and  the  solution  retaining  only  i  =  1  would 
provide  poor  resolution. 

We  evaluate  resolution  of  the  inferred  profile  by  the  general¬ 
ized  Dirichlet  kernel  method  described  in  Sec.  II  and  the  Appendix. 

As  for  the  variance,  the  achievable  height  resolution  depends  on  the 
number  of  spectral  components  retained.  Figure  4  shows  the  modulus 
of  the  averaging  function  A  at  height  44  km  formed  by  truncating  the 
Dirichlet  representation  | Eq .  (25) j  at  various  values  of  k.  The  best 
resolution — i.e.,  the  narrowest  averaging  function — is  achieved  by 
retaining  all  seven  terms.  Unfortunately,  as  discussed  above,  the 
variance  increases  with  the  number  of  terms  retained,  and  a  compromise 
is  needed. 

Figure  5  illustrates  the  tradeoff  between  resolution  and  noise- 
induced  variance  in  the  conductivity  height  profile  computed  for  the 
strong  SPE.  The  resolut ion — def ined  as  the  full  width  at  half-peak 

: k 

This  behavior  occurs  because  higher-order  terms  tend  to  be 
oscillatory  and,  therefore,  exhibit  a  larger  variance. 


Fig.  4--Magnitude  of  the  averaging  function  A  for  the  strong 
SPE  example  (height  =  44  km) 
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of  A — decreases  as  the  number  k  of  retained  terms  is  increased ,  but 

at  the  expense  of  a  greater  variance.  The  vertical  and  horizontal 

1/2 

bars  in  Fig.  3  show  the  resolution  and  [Var  (a)]  ,  respectively 

for  k  =  5,  which  we  chose  as  a  good  compromise.  I f  we  had  retained 
the  maximum  number  of  terms  (k  =  n  =  7),  the  resolution  would  improve 
from  5.1  to  3.5  km,  but  [Var  (a)]  would  increase  by  nearly  two 
orders  of  magnitude.  The  optimum  number  of  terms  depends  on  the  num¬ 
ber  of  frequencies  at  which  reflection  data  are  available,  and  the 
amount  of  measurement  noise.  All  terms  (k  =  n)  should  be  retained 
if  noise  is  negligible,  but  only  a  few  retained  if  the  data  are  noisy. 
The  optimum  value  for  k  must  be  selected  on  a  case-by-case  basis, 
using  tradeoff  graphs  such  as  Fig.  5. 

Physical  Interpretation 

We  next  show  physically  why  the  calculated  profile  shown  in 
Fig.  3  agrees  with  the  true  profile  only  at  heights  between  about  40 
and  50  km.  In  addition,  we  give  diagnostic  methods  for  estimating 
the  height  range  over  which  a  set  of  reflection  data  can  yield  accu¬ 
rate  conductivity  profiles. 

Ground-level  reflection  coefficients  contain  information  about 
only  those  heights  from  which  significant  energy  is  returned.  To 
estimate  the  heights  at  which  important  interactions  between  the  in¬ 
cident  field  and  the  ionosphere  occur,  we  use  the  following  approxi¬ 
mate  expression  (Field  and  Lewinstein  [1978])  for  the  TM  reflection 
coefficient : 


R(0) 


/  L  2ob  S ] 

*'o 


[  —  i 2o>4> ( z )  / c]  , 


(30) 


where 


z 


dz'  q(z’)  , 


(31) 
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2 

and  Q  =  q/n  .  Equation  (30)  is  the  second-order  WKB  approximation, 
and  is  strictly  valid  only  for  weak  reflections.  Nonetheless,  we 
will  show  that  it  provides  useful  insight  into  the  reflection  process. 
The  integrand  %  'of  Eq.  (30)  is  the  contribution  of  a  thin  stratum 
at  height  z  to  the  ground-level  reflection  coefficient.  The  quantity 
in  brackets  represents  conversion  of  upgoing  to  downgoing  waves  in 
the  stratum;  the  exponential  term  represents  the  round-trip  transmis¬ 
sion  loss  suffered  in  propagating  to  and  from  height  z. 

If  we  set 


n  (z)  =  1  -  ia(z)  , 


(32) 


the  reflection-per-unit  height  %  becomes 


#(z)  =  “  - - - — -  exp  [-i2ux| >(z)/c]  (  ,  (33) 

(  (1  -  ia) (C2  -  ia)  j 


which  has  magnitude 


l£(z) 


1  da 

1 

2  2* 
4C  S 

1/2  j 

)  dz  /  4 

,  2  \  1  /  2 

1  ~  ^  2 

( 

(  \c 

+  a  J 

L  1  +  a  J 

) 

exp  [2u)Im  <f>(z)/c] 


(34) 


where 


1.  *U)  -  -  J-  f  ^(c4  +  a2)1/2  -  c2  d,'  . 

*^2  JQ 


(35) 


figure  6  plots  j^|  for  the  exponential  conductivity  profile 
given  by  Eq .  (26)  at  the  frequencies  used  in  the  inversion,  and  Table 
3  summarizes  information  from  the  figure.  The  quantity  zD  denotes 
the  height  at  which  \%\  attains  its  maximum  value,  and  z^  and  z^  de¬ 
note  heights  below  and  above  z^ ,  respectively,  where  |jj;|  is  one-half 
its  peak  value. 


trong  SPE  example 
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Table  3 

REFLECTION  REGIONS  FOR  STRONG  SPE  EXAMPLE 


Frequency 

(kHz) 

ZR 

(km) 

ZB 

(km) 

ZT 

(km) 

Width  (zT  -  zfi) 
(km) 

4 

44.1 

38.8 

50.0 

11.2 

6 

44.1 

39.6 

49.6 

10.0 

9 

44.4 

40.3 

49.3 

9.0 

14 

44.8 

41.1 

49.1 

8.0 

21 

45.5 

41.7 

49.0 

7.2 

30 

45.8 

42.2 

49.0 

6.7 

48 

46.2 

42.7 

49.0 

6.3 

The  maximum  contribution  to  the  ground-level  reflection  co¬ 
efficient  comes  from  heights  near  z  ,  whereas  relatively  little  con- 

R 

tribution  comes  from  heights  below  z^  or  above  z„.  Therefore,  we 
expect  the  ground-level  reflection  coefficient  to  provide  the  best 
information  on  ionospheric  conductivity  in  the  height  range  bounded 
by  and  z and  only  poor  information  outside  of  this  range. 

Table  3  shows  that  z  is  a  slowly  increasing  function  of  fre- 

K 

quency,  but  remains  near  the  center  of  the  40  to  50  km  height  range 
over  the  entire  band  considered.  The  lower  height  z  increases  from 

D 

39  to  43  km  as  the  frequency  is  raised  from  4  to  48  kHz,  but  the  upper 
height  z^,  remains  between  49  and  50  km.  From  Table  3  and  Fig.  6  we 
conclude  that,  for  the  frequencies  used,  the  reflection  coefficients 
give  little  information  about  heights  below  40  or  above  50  km;  and, 
even  within  the  40  to  50  km  height  range,  the  information  is  concen¬ 
trated  between  about  44  and  50  km.  This  conclusion  is  consistent 
with  Fig.  3,  which  shows  excellent  agreement  between  the  calculated 
and  true  profiles  between  44  and  49  km;  fair  agreement  between  40  and 
44  km;  and  a  total  breakdown  below  40  or  above  50  km. 

The  insensitivity  of  z T  to  frequency  indicates  that  the  height 
range  sampled  by  VLF/LF  sounding  could  not  be  expanded  substantirlly 
by  using  higher  frequencies;  and  the  use  of  frequencies  much  lower 
than  4  kHz  is  probably  not  practical.  Thus,  for  the  strong  SPE  ex¬ 
ample,  VLF/LF  sounding  would  appear  incapable  of  determining  the 
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state  of  the  ionosphere  below  40  or  above  50  km.  Finally,  the  fact 
that  z^  is  much  lower  than  70  km  provides  a  posteriori  justification 
of  our  neglect  of  geomagnetic  effects. 


Spacing  of  Frequencies 


We  have  found  that  the  best  results  are  obtained  by  choosing  an 
unequal  spacing  for  the  frequencies  used.  On  physical  grounds  we 
expect  that  the  optimum  choice  of  frequencies  would  correspond  to  a 
uniform  spacing  of  the  heights  of  peak  reflection  z  .  We  therefore 
need  a  mapping  between  reflection  height  and  frequency.  Unfortunately, 
Eq.  (34)  is  inconvenient  for  this  purpose.  However,  Field  and  Engel 
[1965]  showed  that — subject  to  certain  restrictions — the  major  reflec¬ 
tions  take  place  within  a  few  kilometers  of  the  altitude  where 

a  =  Jl  C2  .  (36) 

Specifically,  by  examining  the  second-order  correction  to  the  WKB 
approximation,  they  showed  that  substantial  reflection  can  occur  only 
at  altitudes  where 


2  2 

where  we  have  used  Eq .  (26).  Provided  that  q  <<  n  ,  the  peak  value 
of  the  left  side  of  Eq.  (37)  is  easily  shown  to  occur  at  the  height 
where  a  is  given  by  Eq.  (36).  By  using  Eqs.  (4),  (25),  (32),  and 
(37),  we  find 


r  .  2lTfen"l 

ZR  *  g  \f*  c  ~aW  J  ' 


Table  4  shows  the  values  of  z  given  by  Eq.  (38)  for  the  strong  SPE 
example.  Note  that  the  logarithmically-spaced  frequencies  correspond 
to  evenly-spaced  reflection  heights.  The  differences  between  the 
approximate  z  in  Table  4  and  the  more  accurate  ones  in  Table  3  do 
not  appear  to  strongly  affect  the  optimum  choice  of  frequencies. 


^3=v;‘ 
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Table  4 

APPROXIMATE  PEAK  REFLECTION  HEIGHTS 
FOR  STRONG  SPE  EXAMPLE 


Frequency 

(kHz) 

zR  [Eq.  (38)] 
(km) 

4 

40.7 

6 

41.8 

9 

42.8 

14 

43.9 

21 

45.0 

30 

45.9 

48 

47.1 

In  practice,  we  selected  the  frequency  spacing  for  the  strong 
SPE  model  by  (1)  using  arbitrarily  spaced  frequencies  to  construct 
a  preliminary  profile,  (2)  making  an  exponential  fit  to  this  prelim¬ 
inary  profile,  and  (3)  using  Eq .  (38)  to  select  logarithmically  spaced 

frequencies  corresponding  to  evenly  spaced  z  . 

R 

Height-Gain  Function 

It  is  desirable  to  have  a  method  of  assessing  profiles  as  they 
are  being  constructed.  We  have  found  the  height-gain  function  H(z) 
to  be  useful  in  this  regard.  Here  we  normalize  H  so  that — below  the 
ionosphere — it  equals  the  upgoing  wave,  which  is  normalized  to  unity 
at  the  ground.  The  resulting  expression  is  (Field,  Lewinstein,  and 
Dore  [1976]) 


H 


1  +  R(0) 
1  +  R(z) 


dz’ 


n2(z’)  1 
A(z') 


(39) 


where  R(z)  is  the  reflection  coefficient  and  A(z)  is  the  admittance 


A(z) 


1  +  R(z)  1. 

1  -  R(z)  C  ' 
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Figure  7a  shows  H(z)  for  frequencies  of  4  and  48  kHz  and  the 
exponential  fit  to  the  strong  SPK  |  Eq .  (26) | .  Comparison  with  Fig.  6 
shows  that  most  reflections  occur  in  the  height  range  where  0.5  ^H(z) 

^  1.  At  lower  altitudes  the  ionosphere  is  too  rarefied  to  reflect  the 
wave  strongly;  at  higher  altitudes  the  signal  is  so  weak  that  the  re¬ 
turn  is  small  even  if  the  reflectivity  is  high.  The  shaded  region 
thus  indicates  where  both  the  signal  .inJ  the  reflectivity  are  strong. 
Figure  7b  shows  H  for  the  calculated  profile,  which  we  have  assumed 
retains  its  49-km  value  at  greater  altitudes.  Again,  we  see  that 
0.5  <  H  <_  1  provides  a  measure  of  the  height  range  over  which  the 
calculated  profile  is  valid;  viz.,  between  about  40  and  50  km.  Com¬ 
puting  H  during  profile  construction  gives  a  criterion  for  terminating 
the  iterations.  When  the  major  contribution  for  an  iteration  falls 
outside  the  height  range  defined  by  0.5  <,  H  <,  1,  we  have  done  as  well 
as  possible  and  accept  that  profile  as  the  final  one. 

MODEL  C-LAYER 

The  above  example  pertains  to  conditions  where  the  conductivity 
increases  monotonicallv  with  height.  We  now  consider  a  qualitatively 
different  situation,  where  the  lower  ionosphere  exhibits  a  well-defined 
layer.  Specifically,  we  used  the  model  given  in  Fig.  8  (dashed  lines), 
which  was  shown  by  Rasmussen,  Kossey,  and  Lewis  | 1980]  to  give  close 
agreement  between  measured  and  calculated  ionosound  returns.  This 
profile  consists  of  a  sharply  hounded  layer  of  conductivity  1.8  *  10 
mhos/m  extending  from  63  to  69  km. 

The  magnitudes  of  the  reflection  coefficients  were  measured  at 
an  incidence  angle  of  about  64  deg  by  P.asmussen,  Kossey,  and  Lewis,  and 

k 

were  input  to  our  algorithm  at  10  evenly  spaced  frequencies  between 
5  and  50  kHz.  Since  no  phase  measurements  are  available,  we  had  to 
manufacture  phase  "data"  by  integrating  Eq.  (2)  to  find  Arg  R(0)  for 
the  "true"  profile  in  Fig.  8,  We  used  the  weighting  function  given 
by  Eq.  (26)  with  y  *  0.5  km  A  and  b  =  50  km;  and,  as  before,  we  used 
a  null  (o°  =  0)  initial  profile  and  a  three-point  smoothing  filter. 

* 

The  logarithmic  spacing  discussed  above  is  not  suitable  for  the 
type  of  profile  being  considered  here. 


H  (z)  -•  H  (2) 


Fig.  7--Normal ized  height-gain  function  for  (a)  true 
and  (b)  calculated  profiles 


Height  (km) 


Fig.  8--Calculated  and  true  profiles  after  20  iterations  for  the  C-layer 
example.  Vertical  bars  show  uncertainty  caused  by  hypothetical  noise 
in  the  reflection  coefficients.  Horizontal  bars  show  uncertainty 
caused  by  incompleteness  of  the  reflection  data. 
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Figure  8  shows  the  profile  calculated  from  the  above  reflection 
data  after  20  iterations.  The  calculated  and  true  profiles  agree  well, 
particularly  with  regard  to  the  peak  conductivity  and  location  of  the 
layer.  The  failure  of  the  inverse  solution  to  reproduce  the  sharp 
boundaries  is  caused  by  the  limited  height  resolution  achievable  with 
incomplete  data. 

To  examine  the  effects  of  measurement  noise  on  the  variance  of 
the  computed  profile,  we  use  the  uncertainties  in  reflection  coef¬ 
ficients  given  in  Table  3.  As  before,  these  values  are  representa¬ 
tive  of  actual  measurement  uncertainties.  The  tradeoff  between  vari¬ 
ance  and  height  resolution  at  a  height  of  66  km  is  shown  in  Fig.  9  for 
the  calculated  profile  after  20  iterations.  The  averaging  functions 
for  k  <  4  did  not  display  well-defined  peaks,  so  their  widths  could 
not  be  determined. 


Table  5 

ASSUMED  UNCERTAINTIES  IN  REFLECTION 
COEFFICIENTS  FOR  THE 
C-LAYER  EXAMPLE 


Frequency 

(kHz) 

|S./R.(0)| 

5 

0.3 

10 

0.25 

15 

0.18 

20 

0.15 

25 

0.1 

30 

0.1 

35 

0.1 

40 

0.1 

45 

0.15 

50 

0.15 

The  vertical  and  horizontal  bars  on  Fig.  8  correspond  to  the 

1/2 

uncertainties  caused  by  noise  ([Var  (c)]  )  and  incomplete  data 

(height  resolution)  with  k  =  8.  Note  that  the  agreement  between  the 
calculateu  and  true  profiles  is  much  closer  than  indicated  by  the 
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horizontal  bars.  This  agreement  suggests  that  the  height  resolution 
is,  in  fact,  far  better  than  the  width  of  the  averaging  function  A. 
Shellman  [1973]  also  concluded  that  this  theoretical  estimate  of 
height  resolution  is,  in  many  cases,  much  too  pessimistic.  In  prac¬ 
tice,  the  achievable  resolution  is  probably  best  determined  by  exer¬ 
cising  our  inversion  algorithm  for  a  number  of  generic  profiles,  such 
as  the  two  examples  given  above. 
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IV.  CONCLUSIONS 


Our  extension  of  the  Backus-Gilbert  theory  appears  well  suited 
for  inferring  ionospheric  properties  from  ground-level  VLF/LF  iono- 
sound  data.  It  permits  height  profiles  of  ionospheric  parameters  to 
be  determined  without  resorting  to  trial  and  error  and  explicitly 
gives  the  relationship  between  uncertainties  caused  by  noisy  and  in¬ 
complete  data.  The  analyst  can  therefore  choose  the  best  compromise 
between  variance  and  height  resolution  of  calculated  profiles. 

We  consider  only  altitudes  below  about  70  km,  where  the  propa¬ 
gation  can  be  assumed  isotropic.  This  altitude  restriction  limits 
application  to  certain  ambient  daytime  or  disturbed  conditions.  More¬ 
over,  only  the  profile  for  total  conductivity  can  be  found  from  the 
reflection  data.  Determination  of  profiles  for  the  density  and  col¬ 
lision  frequency  of  each  charged  constituent  requires  auxiliary  in¬ 
formation.  The  method  can  be  extended  to  include  geomagnetic  effects. 

The  versatility  of  our  approach  is  evidenced  by  its  satisfactory 
performance  for  two  qualitatively  different  examples — a  strong  SPE 
where  the  conductivity  increases  quasi-exponentially  with  altitude 
and  a  model  C-layer  exhibiting  a  well-defined  conductivity  peak.  In 
each  example,  the  iterative  solution  converges  toward  the  true  pro¬ 
file  despite  our  use  of  initial  estimates  that  contain  no  information 
about  ionospheric  structure.  The  nonlinear  relationship  between  re¬ 
flection  coefficients  and  ionospheric  conductivity  apparently  causes 
no  difficulties  with  uniqueness. 

We  found  the  height  resolution  of  the  calculated  profiles  for 
the  two  examples  to  be  much  better  than  predicted  by  theory.  Other 
investigators  also  have  reported  that  the  theoretical  limit  on  height 
resolution  is  often  too  pessimistic.  The  resolution  achievable  in 
practice  is  best  determined  empirically  by  inverting  computer-generated 
reflection  coefficients  for  model  conductivity  profiles  that  corre¬ 
spond  to  generic  ionospheric  conditions. 

Ground-level  reflection  data  contain  information  about  only  those 
heights  from  which  significant  energy  is  returned.  Our  calculated 
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profiles  therefore  agree  with  the  true  profiles  only  over  a  limited 
height  range,  which  corresponds  to  the  region  where  the  most  signifi¬ 
cant  reflections  occur.  This  reflection-height  region  provides  a 
physical  basis  for  judging  where  (a)  the  calculated  profiles  are 
reliable,  or  (b)  the  interaction  between  the  ionosound  signal  and 
the  ionosphere  is  too  weak  to  permit  a  valid  inverse  solution. 


1 
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Appendix 


THE  SPECTRAL  EXPANSION  METHOD 


Although  the  inverse  solution  and  its  interpretation  can  be 

carried  out  entirely  in  terms  of  Frechet  kernels  G  (z)  as  discussed 

i 

in  Sec.  II,  significant  advantages  result  from  using  the  set  of  basis 
functions  formed  by  orthogonalizing  the  set  {(J^(z)}.  These  advan¬ 
tages  include  (1)  improved  numerical  stability  in  the  iterative  inverse 
solution,  (2)  the  ability  to  construct  resolution  functions  that  are 
direct  generalizations  of  Dirichlet  kernels  for  Fourier  series,  and 
(3)  explicit  presentation  of  the  tradeoff  between  achieving  high 
resolution  and  minimum  noise-induced  errors.  Our  analysis  has  been 
guided  by  Parker's  [1977]  review,  which  was  based  in  turn  on  work  of 
Gilbert  [1971],  Jackson  [1972],  Wiggins  [1972],  and  Jupp  and  Vozoff 
[1975]. 

The  first-order  variation  in  the  reflection  coefficient  6R^ (0) 
at  frequency  f^,  due  to  a  small  change  in  6o,  is  given  by  Eq.  (11), 
from  which  it  follows  that 


6R.(0) 


S. 

l 


6o(z)  , 


(A .  1 ) 


where  represents  the  measurement  uncertainty  in  R^(0).  Following 
Oldenburg  [1978]  and  Parker  [1977]  we  slightly  generalize  the  matrix 
£  given  by  Eq.  (16)  by  including  the  weighting  function  w(z)  and  the 
error  estimates  S^: 


r 


ij 


dz 


Gi(Z>  G1(z) 


[v(z)]‘ 


(A. 2) 


The  use  of  w(z)  allows  us  to  emphasize  height  regions  where  significant 
conductivity  is  anticipated,  and  suppress  the  nonphysical  build-up  of 
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conductivity  at  low  heights.  By  construction  r  is  a  Hermitian, 
positive-definite  matrix  and  can  therefore  be  diagonalized  by  a 
unitary  matrix  U  as 


u  r  11  =  A 


(A. 3) 


where  (A)..  =  X.6...  Next,  we  define  the  function  H. (z)  as 
~  ij  J  ij  i 


-1/2  V"'  t  Gi  ^ 

V*> s  v 2*<u ’u  s  "<*>  • 

4  J 


(A. 4) 


The  set  {H^(z)}  is  orthonormal  since 


m 

J  dz  Hi(z)H*(z)  =  X“1/2  ^1/22  (U+) 

o  k,e 


ik 


G.  (z)  G. (z)  „ 

x  |  dz  - V-  [w(z)]  U 

Sk  S* 


•C 

/ 


x-l/2  x-V2^(lJt)ik 
k,«, 


.-1/2  ,-1/2  ,  .  =  - 

Xi  Xj  Xi6ij  6ij 


(A. 5) 


The  orthonormality  allows  the  conductivity  correction  6o(z)  to  be 
expanded  as 


n 

&o(z)  -  atH*(z)w(z)  +  6a(z)  , 
1-1 


(A. 6) 


Var  (tf(z)J 


(A.  9) 
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with 


(a .  a.) 

N  i  r 


k,£ 


(A. 10) 


where  the  (  )  denote  ensemble  averages.  Because  are  the  standard 
deviations  in  the  data  and  the  noise  is  assumed  to  be  uncorrelated, 
the  covariance  matrix  for  the  data  is  simply 


and  the  unitarity  of  U  gives 


(A. 11) 


<“i  *j>  •  «ij  •  «-12> 

Substitution  into  Eq.  (A. 9)  produces 

Var  [o(z)l  =  [v(z)]2££  |Hi(z)|2  .  (A. 13) 

i  1 

Equation  (A. 12)  states  that  the  expansion  coefficients  a^  are 
statistically  independent  and  have  variance  1/A^,  where  A^  are  the 
positive-definite  eigenvalues  of  the  T  matrix.  The  spectral  expan¬ 
sion  [Eq.  (A. 6)1  can  therefore  he  interpreted  as  a  superposition  of 
terms  having  progressively  greater  statistical  uncertainty.  The 
di.iRonulization  process  recombines  the  original  Frechet  kernels  Gi 
to  produce  components  H^,  which  allow  the  conductivity  features  that 
can  most  reliably  be  inferred  from  the  data  to  be  identified.  In 
practice,  a  large  reduction  in  Var  (0)  can  be  achieved  by  suppressing 
a  few  of  the  smallest  eigenvalue  terms  from  Eq.  (A. 6).  Because  eigen¬ 
value  components  with  large  At  tend  to  be  smoother  than  those  with 
small  A^,  this  truncation  removes  the  highly  oscillatory  terms,  which 


are  typically  the  noisiest  in  the  data  set.  This  filtering  capability 
is  not  possible  in  the  undiagonalized  basis  set. 

The  resolution  analysis  treated  in  Sec.  II  can  be  assily  carried 
out  in  the  diagonalized  basis.  Denoting  by  (6cf(z))  the  estimate  ol' 
5o(z)  produced  by  retaining  kin  terms  in  the  spectral  expansion 
lEq.  (A. 6)],  we  find 


<6o(z))  =  w(z)  ^  ^  (2) 
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/m  k 

dz'  6o<2,) 

0  i=i 
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dz’  (So(z’)A(z,  z’)  . 


(A. 14) 


The  averaging  function  A  has  the  same  form  as  the  Fourier  series 
Dirichlet  kernel  discussed  in  Sec.  II.  Its  width  is  a  measure  of 
the  maximum  height  range  over  which  the  conductivity  profile  could — 
in  principle — be  deformed  and  still  be  consistent  with  the  reflection- 
coefficient  data  set.  In  practice,  we  have  found  that  the  height 
resolution  is  better  than  indicated  by  this  conservative  criterion. 

In  general,  retaining  fewer  terms  in  F,q.  (A.  14)  has  the  effect 
of  increasing  the  width  of  the  averaging  function.  Because  this 
process  also  decreases  the  variance  in  the  conductivity  [see  Eq .  (A. 13)] 
there  is  a  tradeoff  between  high  resolution  and  st?t:stical  reliability. 
Our  experience  indicates  that  dropping  the  smallest  eigenvalue  terms 
produces  a  significant  improvement  in  estimated  variance,  while  caus¬ 
ing  only  a  negligible  degradation  in  height  resolution.  This  behavior 
is  illustrated  by  the  examples  in  Sec.  III. 


MISSION 

of 

Rome  Air  Development  Center 

MDC  plans  and  execute*  ties  vouch,  development,  test  and 
selected  acquisition  programs  in  support  o£  Command ,  Contact 
Communications  and  Intelligence  (C3I)  activities.  Technical 
and  engineering  support  within  areas  of  technical  competence 
as  provAded  to  ESI?  Program  Ofoices  (Ft?*)  and  other  ESP 
elements.  The  principal  technical  mission  areas  are 
communications,  electromagnetic  guidance  and  control,  sur¬ 
veillance  £>tj  ground  and  aerospace  objects,  intelligence  data 
collection  and  handling ,  information  system  technology, 
ionospheric  propagation ,  solid  state,  sciences ,  microwave 
physics  and  electronic  reliability,  maintainability  and 
compatibility. 


